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In a number of astrophysical applications one tries to determine the two- 
dimensional or three-dimensional structure of an object from a time series of 
measurements. While most methods used for reconstruction assume that ob- 
ject is static, the data are often acquired over a time interval during which the 
object may change significantly. This problem may be addressed with time- 
dependent reconstruction methods such as Kalman filtering which models the 

£j , temporal evolution of the unknown object as a random walk that may or may 

c/3 ■ 

ci . not have a deterministic component. Time-dependent reconstructions of a hy- 

j> ! drodynamical simulation from its line-integral projections are presented. In these 

^ , examples standard reconstructions based on the static assumption are poor while 

the Kalman based reconstructions are of good quality. Implications for various 
astrophysical applications, including tomography of the solar corona and radio 
aperture synthesis, are discussed. 



1. Introduction 

Tomographic reconstruction has a number of astrophysical applications such as de- 
termining the three-dimensional structure of the solar corona and heliosphere with radio 
scintillation, white-light or EUV data (Dunn et al. 2005, Frazin 2000, Frazin & Janzen 2002, 
Frazin et al. 2005, Hayashi et al. 2003), Doppler tomography of accreting flows (e.g., Marsh 
2005), air-column tomography in adaptive optics (e.g., Tyler 1994, Tokovinin et al. 2001), 
and time-distance helioseismic tomography of near-surface features (e.g., Duval & Gizon 
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2000) . It is quite often the case that the object exhibits significant evolution during the time 
in which the data required for the inversion are collected. For example, white-light tomog- 
raphy of the solar corona requires 180° of rotation, which takes about two weeks. During 
some two week time periods the large scale topology of the structure may not change a lot, 
but during others it can change dramatically, which affects the quality of the reconstructions 
(Butala et al. 2005). Using two or more spacecraft has been shown to improve the situa- 
tion, but the problem demands tomographic inversion methods which explicitly account for 
temporal variations (Frazin & Kamalabadi 2005). 

Determining the time-dependent structure of an object with tomographic methods is 
an example of what is sometimes called dynamic estimation in the signal processing and 
statistical literature. Dynamic estimation also has important applications in radio astronomy 
imaging. Standard approaches to radio-interferometric image formation (Thompson et al. 

2001) typically assume the source brightness distribution under reconstruction has negligible 
time-variability over the course of the observing period. There are astrophysical sources for 
which this standard assumption in Fourier synthesis imaging cannot be made. Studies of 
highly time- variable solar phenomena (Bastian 1989), comets colliding with planets (de Pater 
& Brecht 2001), and relativistic jet sources ( "microquasars" ) in the Galaxy (Vermeulen 1993; 
Mioduszewski 2001) may show proper motions or flux density variability within the span of 
a given observing run. Dynamic estimation offers the possibility of improved image fidelity 
in radio-interferometric imaging of objects of this type. Historically, several approaches have 
been taken to mitigate this problem, including sub-dividing the measured visibility data into 
snapshot time intervals, over which the standard assumption of a constant source brightness 
distribution has greater validity. These data segments can then be imaged separately, but 
at the cost of reduced image fidelity as a result of sparser u-v plane coverage and lower 
sensitivity in each individual interval. Furthermore, this approach does not take advantage 
of the fact that an image at one measurement time is closely related to the image at the 
next measurement time. 

Dynamic estimation is closely related to data assimilation. Data assimilation incorpo- 
rates physics-based models for the temporal evolution and is commonplace in the atmospheric 
science, weather prediction, oceanographic and other communities (e.g., Ghil 1989, Seppanen 
et al. 2001, Bertino et al. 2002, Buehner & Malanotte-Rizzoli 2003). The problem addressed 
here is that of dynamic estimation when no physics-based model for the temporal evolution 
is available. 

On the surface, it may seem that the task of determining a three-dimensional object (2 
spatial dimensions plus time) from two dimensions of data (1 spatial dimension plus rotation 
angle) in the case of interferometry, or a four-dimensional object from three dimensions of 



- 3- 



data in the case of tomography, would be impossible to do in a useful way. Even the static 
reconstruction problem is often somewhat underdetermined or ill-conditioned, let alone the 
time-dependent problem. It is helpful to recall that the reason the static reconstruction 
problems are solvable is because real objects are smooth or can be approximated reasonably 
with smooth objects (possibly allowing for edges and ridges). Similarly, a dynamic object 
tends to change in a smooth way. Thus, many real dynamic objects can be usefully approx- 
imated by smooth objects that vary smoothly in time. Smoothly evolving smooth objects 
exist in a much smaller subspace than is available to arbitrary objects that need not ful- 
fill any constraints. The Kalman filter approach presented below effectively finds the most 
smoothly evolving smooth object that matches the data. 

Dynamic tomography without physics-based models has received some attention in the 
literature and various methods have been proposed and demonstrated. Some methods are 
based on Kalman-type filters (e.g., Vauhkonen et al. 2001, Kolehmainen et al. 2003) while 
others are not (Schmitt & Louis 2002, Schmitt et al. 2002, Zhang et al. 2005). There is also a 
large literature on handling patient movement and quasi-periodic cardiac motion in medical 
tomography that can be found in the IEEE Transactions on Medical Imaging. The results 
given this letter are based on the Kalman filter formulation described in the next section. 



2. Time-Dependent Image Reconstruction with the Kalman Filter 

Kalman filtering has been used for time-dependent estimation problems (such as tracking 
satellites) since its introduction in 1960 (e.g., Kalman 1960; Anderson & Moore 1979; Tapley 
et al. 2004). The Kalman filter is based on the so-called state-space formulation, which 
consists of a time-update equation: 

x t+1 = U t x t + g t , (1) 

and a measurement equation: 

y t = A t x t + n t , (2) 

where t is a time index, and the vector x t is a discrete representation of the unknown 
object (e.g., the electron density of the solar corona) at time t. The matrix Ut is called the 
update operator, A t is the measurement operator for all measurements made at (or suitably 
close to) time t, and y t is the vector of measurements taken at time t (e.g., a polarized 
brightness image from a coronagraph) . The vector x t is also called the state of the system 
at time t. The Gaussian zero-mean vector n t represents noise in the measurements, and 
the Gaussian zero-mean vector g t , called the state noise, represents a noise process that 
drives non-deterministic changes in the state from one time to the next. It is the state 
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noise that gives the model its random walk character and its smooth time evolution. In 
data assimilation the update operator U t contains knowledge of the physics of the system 
(e.g., a hydrodynamical model) and is responsible for the deterministic character of the 
system behavior. When no good physics model is available U t is usually just the identity 
operator, but it can also be used to perform simple deterministic tasks, such as applying 
a differential rotation curve in the case of solar tomography. The Kalman formulation also 
requires knowledge of the covariance matrices g k gf, riknf, and g^nf, where k and I are time 
indices (the T superscript is the transpose operator and the overline represents averaging over 
the statistical ensemble). Kalman filtering allows one to make an estimate, x t \ t > : of the state 
at any time step t given data up to and including time step t', where t' may be greater than 
("smoothing"), less than ("predicting"), or equal to ("filtering") t. 

When solving underdetermined or poorly conditioned linear systems of equations it 
is common practise to regularize the problem (Karl 2000, Frazin & Janzen 2002) to avoid 
spurious, high-frequency oscillations in the solution. Regularization creates smooth solutions. 
Most treatments of the Kalman filter do not incorporate regularization, but our simulations 
were much improved with regularization. Following Baroudi et al. (1998), we regularized the 
solution by augmenting the measurement equation. Instead of measurement equation (2), 
we used the following augmented system: 



Vt 




At 
R 



where is a vector of zeros with the appropriate length and R is a smoothing matrix, which 
was taken to a finite different approximation of the gradient operator. The vector h is a 
noise process. It is assumed that n t h T = 0, and the covariance matrices of n t and h must 
be specified. They will be denoted as C n = n t nf and Ch = hh T , respectively. Rigorously, 
the positive definite matrix C n should include pixel-to-pixel correlations in observational 
errors, but the most important aspect is that the diagonal elements are set to the proper 
noise level. The matrix Ch may take any positive definite form. For the examples below 
we used Ch = o" 2 l, where 1 is the identity matrix, and a is a regularization parameter. To 
understand the effect of this regularization parameter, consider the weighted least-squares 
solution of equation (3) under the assumption Gaussian noise. The solution is given by the 
minimum of the following function: $(cc t ) = (y t — A t x t ) T C n 1 (y t — A t x t ) + (l/a 2 )xj R T Rx t 
(Moon & Sterling 2000, Frazin & Kamalabadi 2005). This is the standard form for a static, 
regularized least-squares solution given only the data y t . When a is small the solution will 
be highly smoothed and the data misfit vector (y t — A t x t ) will be large. Thus, a serves as 
a regularization parameter which may be chosen via methods such as cross-validation (Karl 
2000, Frazin & Janzen 2002, and references therein). 
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In addition, we were able to improve the solution by adopting a non-diagonal form 
of the state noise covariance matrix C g = g t gf, which must be done with care so as to 
ensure that it is positive definite. In the simulations presented below, we set C g so that 
the correlation in the state noise between any two pixels decreased exponentially with the 
square of the distance between them. The size of the exponential decrease was characterized 
by a correlation length I. 

The formulae for the estimates x t ui can be found in standard texts such as Anderson 
& Moore (1979), Kailath et al. (2000) and Tapley et al. (2004). For numerical stability, 
all results given below were calculated with the so-called "square-root" or "array" form of 
the Kalman smoother. All reconstructions given in the next sections are estimates based on 
all of the data, i.e., the dynamic reconstructions shown below are elements of the sequence 
x t \N, where t is the time index of the image and N is the total number of observation times. 
These are the so-called "smoothed" estimates in the literature because N > t (except for 
the last estimate x N \ N ), which has little to with smoothing in the regularization sense. 

3. Simulation Results 

In order to demonstrate the principles given above we made tomographic reconstructions 
of a highly time-variable two-dimensional (2D) object from its line-of-sight projections. Note 
that this problem considered here is equivalent to aperture synthesis of a circumpolar source 
with a linear array (Bracewell & Riddle 1967). The time- variable object was given by a 
2D MHD simulation of a magnetized molecular cloud collapse by C.F. Gammie and C. 
Ditsworth. 1 This is a highly dynamic simulation with iV = 128 time-steps or movie frames. 
At the beginning of the simulation the cloud is uniform with sinusoidal density and magnetic 
field perturbations superimposed. As the simulation progresses, the cloud collapses into 
several strands which begin to orbit around each other as they each undergo further collapse 
to form protostars by the final time-step. 

The simulated data were given by line-of-sight integrals of the MHD movie frames. 
Between each MHD time-step, the view-angle was increased by 5.6°, so that every 64 frames 
the view angle cycles through 360°, making 2 complete rotations during the simulation 
period. Since data taken at an angle 9 and 9 + 180° are redundant, 32 projections (one per 
time-step) were needed to make a static reconstruction with full angular coverage. The first 
half of the movie shows the object undergoing relatively slow changes, while the second half 



1 The 2D cloud collapse simulation can be found on the World Wide Web at: 
http: / / ddr.astro.uiuc.edu/ddr/twod. 



- 6- 



contains most of the rapid dynamics. White, Gaussian noise was added to the projections 
such that y had a signal-to-noise ratio (SNR) of 5 x 10 3 , where SNR = y/a n . 

Figure 1 shows the results of the simulations. Each row represents a different time 
in the MHD movie. The first column shows the density of the MHD movie at the time 
step m (indicated in the figure), the second column shows the Kalman reconstruction x m \N, 
and the third column shows the static reconstruction based on 32 frames made with data 
centered on that time step. The state noise correlation length /, state noise amplitude 
yjL(C g ), and regularization parameter a where all chosen "by eye" to give the best looking 
reconstructions. The "optimal" parameters were: /=3.6 pixels, ^L(C g ) = 5 x 10 -3 af, 
and cr = 0.02 a n L(R)/L(A ), where L(-) indicates that the largest singular value of the 
matrix argument is to be taken. More rigorous methods such as cross-validation may also 
be used to determine reconstruction parameters (Karl 2000, Frazin & Janzen 2002, and 
references therein). For the static reconstructions, the regularization parameter was chosen 
the minimize the difference between the state estimate and the true state, averaged over the 
32 frames used to form each static reconstructions. 

The results presented in Figure 1 show that the Kalman reconstructions are much more 
faithful to the original than the static reconstructions. More quantitatively, Figure 2 provides 
a comparison between the Kalman and static reconstruction error, e m . In the Kalman case, 
e m is defined by: e m = \\x m — cc m |7v||, with an analogous definition for the static error. At 
times greater than m = 66 the MHD simulation begins to exhibit extremely rapid temporal 
changes, and both the Kalman and static reconstructions match the simulation frames poorly. 
However, the results indicate that in the regime of low to moderate temporal variability, the 
Kalman reconstructions are far superior. Quantifying the effects of temporal variability is 
a subject of ongoing research. Note that the far superior qualitative morphology of the 
Kalman over the static reconstructions evident in Figure 1 does not necessary result in a 
far superior error in Figure 2. The quantification of reconstruction quality is also a topic of 
continuing study. 

The left plot of Figure 3 provides a comparison between the p'th power of the Kalman 
error (with p = 2.5) and a measure of the time variability in the MHD simulation, d m , defined 
by the centered difference: d m = \ ||cc m +i — cc m _i||. The Kalman error has been raised to the 
p'th power to aid in the empirical agreement between the two quantities. Both (e m ) p and d m 
have been normalized by dividing by their respective maximum for the purpose of displaying 
them simultaneously in the left plot of Figure 3. The right plot in Figure 3 shows log(d m ) 
vs. log(e m ) and the line of best fit. From Figure 3, it is clear that (e m ) p and d m are highly 
correlated, thus showing the relationship between temporal variability and reconstruction 
error. 
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4. Conclusion 

This paper has demonstrated promising results for the Kalman solution of dynamic 
inverse imaging problems in which the object is evolving in time as the data are being 
collected. It was shown that the Kalman solutions were superior to static reconstructions 
(based on sliding time windows) except for time periods of extremely rapid variability, in 
which case both methods were equally inadequate. The Kalman reconstructions presented 
here were implemented with two features to improve the solutions: 

1. The solutions was regularized by replacing equation (2) with (3), which forces smooth 
solutions. 

2. The random walk was given a correlation length with a non-diagonal form of the state 
noise covariance matrix C g . 

There is much to be done along this line of research. Some of the significant issues that 
need to be addressed are: 

• Enforcing positivity constraints: For many astrophysical applications the unknown 
quantity to be imaged must be positive. One approach for this is suggested by Simon 
& Chia (2002). 

• Reducing computational complexity: In order to solve large 2D and 3D problems, the 
computational complexity must be reduced. Much work along these lines has already 
been done (e.g., Cane et al. 1996, Asif & Moura 1999, Treebushny & Madsen 2003, 
Khellah et al. 2005, and references therein) and it is an ongoing area of investigation. 

• Assimilating physics-based models: For solar work, using an MHD model for the update 
operator in equation (1) will greatly increase the ability of the filter to track temporal 
changes. For microquasar imaging in radio interferometry a kinematic model may be 
suitable. 

The authors thank C.F. Gammie for the MHD simulation. This research was supported 
by NASA Sun-Earth Guest Investigator Program grant NNG04GG32G to the University of 
Illinois. 
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Fig. I. — A comparison between the static and dynamic reconstructions of several frames of 
a 2D MHD simulation of a self-gravitating, magnetized cloud. The above panel is divided 
into three rows, each associated with a time index m of the 2D MHD simulation. All images 
in a given row are displayed using the same color scale as the 2D MHD simulation frame of 
that row. The columns of the above panel, from left to right, contain: frames from the 2D 
MHD simulation, Kalman reconstructions, and static reconstructions. Note the superiority 
of the Kalman over the static reconstructions (in particular, the static reconstructions are too 
smooth). In both cases, the quality of the reconstructions decreases as the time- variability 
in the frames of the 2D MHD simulation increases. 
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Fig. 2. — A comparison between the Kalman and static reconstruction. The Kalman recon- 
structions are superior to the static ones, except above time m = 55 where the object begins 
to change very quickly (in which case the two methods have equal errors). 
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Fig. 3. — Left: The normalized magnitude of the dynamic estimation error to the p'th 
power (with p = 2.5) and the normalized magnitude of the centered difference of the state 
versus the time index. The power of p was chosen to aid in the agreement between the two 
quantities. It is clear that the two quantities show a strong correlation. Right: A log-log 
plot of the centered time difference vs. the dynamic estimation error. 



